Hacer una primera aproximación a los datos disponibles sobre la contaminación en el municipio de Madrid.
Vamos a analizar los datos de NO2 de 2018 siguiendo estos pasos:
Descarga y extracción de los datos
Limpieza de de los datos
Análisis univariante de la serie temporal
Enriquecimiento de los datos con otras fuentes
Análisis multivariante de los datos
Planteamiento de un modelo predictivo
Las librerías que vamos a usar son:
pkgs <- c("lubridate", "data.table", "ggplot2", "prophet", "ggExtra", "leaflet", "magrittr", "dplyr", "leaflet", "gplots", "xts")
lapply(pkgs, library, character.only = TRUE, quietly = T)
Los datos de calidad del aire se han descargado desde el portal de datos abiertos de Madrid y se han almacendo descomprimidos dentro de la carpeta Anio201810 en data_input.
Vamos a unificar la información de los ficheros de la calidad del aire en un único data frame llamado raw_air_quality.
path_files_air <- "data_input/Anio201810/"
file_list <- list.files(path = path_files_air, pattern = ".csv")
raw_air_quality <- data.frame()
for (file in 1:length(file_list)) {
path_file <- paste0(path_files_air, file_list[file])
raw_file <- read.csv(path_file, header = TRUE, sep=";")
raw_air_quality <- rbind(raw_air_quality, raw_file)
}
rm(raw_file)
dim(raw_air_quality)
## [1] 54831 56
Cada fila hace referencia a todas las medidas tomadas hora a hora durante el mismo día en una misma estación.
names(raw_air_quality)
## [1] "PROVINCIA" "MUNICIPIO" "ESTACION" "MAGNITUD"
## [5] "PUNTO_MUESTREO" "ANO" "MES" "DIA"
## [9] "H01" "V01" "H02" "V02"
## [13] "H03" "V03" "H04" "V04"
## [17] "H05" "V05" "H06" "V06"
## [21] "H07" "V07" "H08" "V08"
## [25] "H09" "V09" "H10" "V10"
## [29] "H11" "V11" "H12" "V12"
## [33] "H13" "V13" "H14" "V14"
## [37] "H15" "V15" "H16" "V16"
## [41] "H17" "V17" "H18" "V18"
## [45] "H19" "V19" "H20" "V20"
## [49] "H21" "V21" "H22" "V22"
## [53] "H23" "V23" "H24" "V24"
Vamos a trabajar sólo con los datos de NO2, cuya magnitud está codificada como 8.
air_quality <- raw_air_quality[raw_air_quality$MAGNITUD == 8, ]
Para simplificar posteriores análisis, añadimos una nueva columna con el código de las estaciones que corresponde con los 8 primeros dígitos del campo PUNTO_MUESTREO. Es la concatenación del los identificadores de provincia más municipio junto con el de la estación.
station <- as.numeric(substr(air_quality$PUNTO_MUESTREO, 1, 8))
Reajustemos el formato para que cada línea contenga sólo datos verificados, por estación, día y hora.
Vamos a usar un data frame, air_quality_by_days, con la información de estación y fecha que utilizaremos como base para desglosar por fila los valores de NO2 y validez, por cada hora del día.
air_quality_by_days <- data.frame(
station, year = air_quality$ANO,
month = ifelse(air_quality$MES < 10, paste0("0", air_quality$MES), air_quality$MES),
day = air_quality$DIA,
stringsAsFactors = FALSE
)
str(air_quality_by_days)
## 'data.frame': 8750 obs. of 4 variables:
## $ station: num 28079004 28079004 28079004 28079004 28079004 ...
## $ year : int 2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
## $ month : chr "04" "04" "04" "04" ...
## $ day : int 1 2 3 4 5 6 7 8 9 10 ...
Nota: Resto una unidad al índice de hora para que la hora esté en el rango 0 a 23.
air_quality_by_hour <- data.frame()
for (hour in 1:24) {
air_quality_by_days$hour <- hour - 1
air_quality_by_days$no2 <- air_quality[ , (2 * hour + 7)]
air_quality_by_days$verif <- air_quality[ , (2 * hour + 8)]
air_quality_by_hour <- rbind(air_quality_by_hour, air_quality_by_days)
}
rm(air_quality_by_days)
Nos quedamos sólo con los datos verificados:
air_quality_by_hour$no2[air_quality_by_hour$verif != "V"] <- NA
Y añadimos una columna que contenga la fecha en formato estándar.
air_quality_by_hour$date <- ymd_hms(paste0(
air_quality_by_hour$year, "-", air_quality_by_hour$month, "-", air_quality_by_hour$day,
" ", gsub(pattern = "H", replacement = "", air_quality_by_hour$hour), ":00:00"
))
str(air_quality_by_hour)
## 'data.frame': 210000 obs. of 8 variables:
## $ station: num 28079004 28079004 28079004 28079004 28079004 ...
## $ year : int 2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
## $ month : chr "04" "04" "04" "04" ...
## $ day : int 1 2 3 4 5 6 7 8 9 10 ...
## $ hour : num 0 0 0 0 0 0 0 0 0 0 ...
## $ no2 : num 21 67 14 8 20 71 32 23 13 7 ...
## $ verif : Factor w/ 2 levels "N","V": 2 2 2 2 2 2 2 2 2 2 ...
## $ date : POSIXct, format: "2018-04-01" "2018-04-02" ...
Primero obtengamos una visión de la tendencia en 2018:
air_quality_by_day <- air_quality_by_hour
air_quality_by_day$date <- as.Date(air_quality_by_day$date, format = "%Y-%m-%d")
dt_no2 <- data.table(air_quality_by_day)[ , .(median_no2 = median(no2, na.rm = T)), by = .(date)]
names(dt_no2) <- c("ds", "y")
ggplot(data = dt_no2, aes(x = ds, y = y)) + geom_line(color = "#1F77B4") +
geom_smooth(method = 'loess', color = "#1F77B4", linetype = 2) +
labs(x = "", y = expression(NO[2]))
Y por días de la semana:
air_quality_by_day$weekday <- weekdays(air_quality_by_day$date)
air_quality_by_day$weekday <- factor(
air_quality_by_day$weekday,
levels=c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")
)
air_quality_by_day <- air_quality_by_day[order(air_quality_by_day$weekday), ]
ggplot(air_quality_by_day, aes(x = weekday, y = no2)) + geom_boxplot(fill = "#5fa2ce") +
labs(x = "", y = expression(NO[2]))
## Warning: Removed 919 rows containing non-finite values (stat_boxplot).
Obsevando la componente semanal y la anual podemos llegar a la conclusión de que el nivel de NO2 baja en días festivos, pero también en los no lectivos.
Como boceto visual de cómo serían los mapas de calor muestro éste de horas frente días del mes (nota: habría que eliminar los datos no verificados y cambiar la leyenda de los ejes para que fuera más fácil de leer).
air_quality <- select(raw_air_quality, -contains("V"))
air_quality_by_day_by_hour <- select(air_quality, 7:ncol(air_quality))
aggdata <- aggregate(
air_quality_by_day_by_hour,
by = list(air_quality_by_day_by_hour$DIA),
FUN = mean, na.rm = TRUE
)
aggdata <- data.matrix(aggdata[ -c(1, 2)])
colheatmap <- c("#1c5998", "#1c73b1", "#3a87b7", "#67add4","#7bc8e2","#cacaca", "#fdab67", "#fd8938","#f06511", "#d74401", "#a33202")
heatmap.2(
aggdata,
dendrogram = "none",
Colv = NA,
Rowv = NA,
density.info = "none",
trace = "none",
col = colheatmap
)
Como era previsible, en las horas nocturnas, con poco tráfico rodado, también baja el nivel de NO2. Quedaría pendiente analizar, desglosándolo por meses, por qué a veces se produce una mayor concentración a última hora del día, sobre todo a final de mes.
¿Y si vemos cómo se comportan los valores máximos y mínimos frente a la mediana?
dt_no2 <- data.table(air_quality_by_day)[ , .(
median_no2 = median(no2, na.rm = T),
min_no2 = min(no2, na.rm = T),
max_no2 = max(no2, na.rm = T) ), by = .(date)]
ggplot(dt_no2, aes(date)) +
geom_line(aes(y = min_no2), color = "#5f9ed1") +
geom_line(aes(y = max_no2), color = "#fbb04e") +
geom_line(aes(y = median_no2), color = "#006ba4") +
labs(x = "", y = expression(NO[2]))
Tanta diferencia entre los valores extremos dentro del mismo día pueden ser debidos a que las mediciones se han tomado en distintas zonas geográficas. Debemos introducir más fuentes de datos para poder seguir investigando.
He probado Prophet, pero en este caso no lo recomiendo pues hace una previsión de NO2 muy alta para enero de 2019, lo que desvirtúa el resto del análisis. Posiblemente se corregiría si incluyéramos más años en el análisis, o bien podríamos hacer un análisis ARIMA.
dt_no2 <- data.table(air_quality_by_day)[ , .(median_no2 = median(no2, na.rm = T)), by = .(date)]
names(dt_no2) <- c("ds", "y")
m <- prophet(dt_no2, daily.seasonality = TRUE, yearly.seasonality = TRUE )
future <- make_future_dataframe(m, periods = 360)
forecast <- predict(m, future)
plot(m, forecast)
prophet_plot_components(m, forecast)
Otros análisis mulltivariantes que podrían ser interesantes:
Vamos a enriquecer estos datos con la geolocalización de las estaciones de medición, así como la medida de la temperatura.
Utilizaremos la información de las estaciones contenida en el fichero madrid_air_quality_stations.csv de la carpeta data_input, extraído de aquí:
stations <- read.csv("data_input/madrid_air_quality_stations.csv", header = TRUE, sep = ",")
str(stations)
## 'data.frame': 24 obs. of 8 variables:
## $ station : int 28079004 28079008 28079011 28079016 28079017 28079018 28079024 28079027 28079035 28079036 ...
## $ area : Factor w/ 5 levels "centro","noreste",..: 1 1 1 2 5 5 3 2 1 4 ...
## $ name : Factor w/ 24 levels "Arturo Soria",..: 17 10 2 1 24 11 5 3 18 14 ...
## $ address : Factor w/ 24 levels "Avd. Betanzos esq. C/ Monforte de Lemos",..: 21 17 6 8 13 15 16 14 22 2 ...
## $ altitude : int 635 670 708 693 604 630 642 621 659 685 ...
## $ type : Factor w/ 3 levels "S","UF","UT": 3 3 3 2 2 2 1 2 2 3 ...
## $ longitude: num -3.71 -3.68 -3.68 -3.64 -3.71 ...
## $ latitude : num 40.4 40.4 40.5 40.4 40.3 ...
df_air <- merge(air_quality_by_hour, stations)
str(df_air)
## 'data.frame': 210000 obs. of 15 variables:
## $ station : num 28079004 28079004 28079004 28079004 28079004 ...
## $ year : int 2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
## $ month : chr "04" "04" "04" "04" ...
## $ day : int 1 2 3 4 5 6 7 8 9 10 ...
## $ hour : num 0 0 0 0 0 0 0 0 0 0 ...
## $ no2 : num 21 67 14 8 20 71 32 23 13 7 ...
## $ verif : Factor w/ 2 levels "N","V": 2 2 2 2 2 2 2 2 2 2 ...
## $ date : POSIXct, format: "2018-04-01 00:00:00" "2018-04-02 00:00:00" ...
## $ area : Factor w/ 5 levels "centro","noreste",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 24 levels "Arturo Soria",..: 17 17 17 17 17 17 17 17 17 17 ...
## $ address : Factor w/ 24 levels "Avd. Betanzos esq. C/ Monforte de Lemos",..: 21 21 21 21 21 21 21 21 21 21 ...
## $ altitude : int 635 635 635 635 635 635 635 635 635 635 ...
## $ type : Factor w/ 3 levels "S","UF","UT": 3 3 3 3 3 3 3 3 3 3 ...
## $ longitude: num -3.71 -3.71 -3.71 -3.71 -3.71 ...
## $ latitude : num 40.4 40.4 40.4 40.4 40.4 ...
Utilizaremos la información de temperaturas medias por días contenidas en el archivo madrid_hourly_temperatures_2018.csv de la carpeta data_input, extraído de aquí.
Llamaremos date_day a la columna con la información del día y date a una nueva columna que incluya también la hora:
temperatures <- read.csv("data_input/madrid_hourly_temperatures_2018.csv", header = TRUE, sep = ",")
temperatures$date_day <- ymd(temperatures$date)
temperatures$date <- ymd_hms(paste(temperatures$date_day, temperatures$hour), truncated = 3)
str(temperatures)
## 'data.frame': 8760 obs. of 4 variables:
## $ date : POSIXct, format: "2018-01-01 00:00:00" "2018-01-01 01:00:00" ...
## $ hour : int 0 1 2 3 4 5 6 7 8 9 ...
## $ temp : num 7.3 7.4 6.8 7.1 5.3 3.8 2.9 2.9 2.3 4.2 ...
## $ date_day: Date, format: "2018-01-01" "2018-01-01" ...
Guardamos toda la información en un único data frame y lo exportamos a un fichero para futuros análisis:
df_air <- merge(df_air, temperatures)
write.csv(df_air, file = "data_output/df_air.csv", row.names = FALSE)
str(df_air)
## 'data.frame': 210000 obs. of 17 variables:
## $ hour : num 0 0 0 0 0 0 0 0 0 0 ...
## $ date : POSIXct, format: "2018-01-01" "2018-01-01" ...
## $ station : num 28079060 28079038 28079008 28079057 28079055 ...
## $ year : int 2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
## $ month : chr "01" "01" "01" "01" ...
## $ day : int 1 1 1 1 1 1 1 1 1 1 ...
## $ no2 : num 8 13 33 5 9 12 5 32 37 8 ...
## $ verif : Factor w/ 2 levels "N","V": 2 2 2 2 2 2 2 2 2 2 ...
## $ area : Factor w/ 5 levels "centro","noreste",..: 2 1 1 2 2 1 2 5 5 4 ...
## $ name : Factor w/ 24 levels "Arturo Soria",..: 21 7 10 20 22 4 12 11 19 9 ...
## $ address : Factor w/ 24 levels "Avd. Betanzos esq. C/ Monforte de Lemos",..: 23 5 17 11 12 1 18 15 24 3 ...
## $ altitude : int 715 698 670 700 618 674 660 630 604 627 ...
## $ type : Factor w/ 3 levels "S","UF","UT": 2 3 3 2 2 3 1 2 3 2 ...
## $ longitude: num -3.69 -3.71 -3.68 -3.66 -3.58 ...
## $ latitude : num 40.5 40.4 40.4 40.5 40.5 ...
## $ temp : num 7.3 7.3 7.3 7.3 7.3 7.3 7.3 7.3 7.3 7.3 ...
## $ date_day : Date, format: "2018-01-01" "2018-01-01" ...
Al analizar la relación entre el nivel de NO2 y la temperatura no se observa que exista correlación:
p <- ggplot(df_air, aes(x = temp, y = no2, color = -temp)) +
geom_point() +
theme(legend.position = "none") +
labs(x = "Temperatura", y = expression(NO[2]))
ggMarginal(p, type = "histogram", fill = "#1F77B4")
## Warning: Removed 919 rows containing missing values (geom_point).
Ni aunque lo desglosemos por cada una de las estaciones:
ggplot(df_air, aes(x = temp, y = no2, color = -temp)) +
geom_point() +
theme(legend.position = "none") +
facet_wrap( ~ name) +
labs(x = "Temperatura", y = expression(NO[2]))
## Warning: Removed 919 rows containing missing values (geom_point).
Lo que sí nos da una pista es que existe mucha variación entre el nivel de NO2 de las distintas estaciones.
Vamos a comparar los valores medianos haciendo tres grupos según al cuartil en el que estén los datos:
station_media <- df_air %>%
group_by(name) %>%
summarise( no2 = median(no2, na.rm = TRUE)) %>%
arrange(desc(no2))
station_media <- as.data.frame(station_media)
station_media$name <- factor(station_media$name,
levels = station_media$name)
station_media <- station_media[order(station_media$name), ]
station_media$level_no2 <- "high"
station_media$level_no2[station_media$no2 < quantile(station_media$no2)[4]] <- "medium"
station_media$level_no2[station_media$no2 < quantile(station_media$no2)[2]] <- "low"
ggplot(station_media, aes(x=name, y = no2, group = factor(level_no2), fill = factor(level_no2))) +
scale_fill_manual(values = c("#fc7d0b", "#1170aa", "#5fa2ce")) +
geom_bar(stat = "identity") +
coord_flip() +
theme_minimal() +
theme(legend.position = "none") +
labs(x = "", y = expression(NO[2]))
Y la relación entre las estaciones y su localización geográfica:
station_media <- merge(station_media, stations)
station_media$color_no2 <- "red"
station_media$color_no2[station_media$no2 < quantile(station_media$no2)[4] ] <- "orange"
station_media$color_no2[station_media$no2 < quantile(station_media$no2)[2] ] <- "green"
icons <- awesomeIcons(
icon = 'ios-close',
iconColor = station_media$color_no2,
library = 'ion',
markerColor = station_media$color_no2
)
leaflet(station_media) %>% addTiles() %>%
addAwesomeMarkers(lng = station_media$longitude,
lat = station_media$latitude,
icon = icons,
label = ~paste(name, "-", as.character(no2)))
Podemos apreciar cómo las zonas más céntricas contienen mayor nivel de NO2. La excepción es la estación que se encuentra situada en el Parque de «El Retiro», que tiene un nivel de 22, menos de la mitad que «Escuelas Aguirre» con 51.
Nota: A la vista de los resultados deberíamos incluir «Castellana» dentro del grupo high.
Variable a predecir: NO2
Variables independientes:
Sería deseable tener variables meteorológicas por cada estación por día e incluso por hora.
Además, tener histórico de más años, nos va a ayudaría a determinar la estacionalidad de los datos.
Primero dividiríamos el conjunto de datos en dos grupos: entrenamiento y testeo.
Para evitar sobreajustes podemos personalizar los parámetros del conjunto de entrenamiento con la función trainControl de la librería caret, por ejemplo, haciendo k-fold cross-validation que se repita r veces
library(caret)
control <- trainControl(method = "repeatedcv", number = k, repeats = r)
Podríamos probar con los algoritmos: randomForest, xgboost, glmnet para ver cuál hace una mejor predicción:
library(randomForest)
library(xgboost)
library(glmnet)
fit.rf <- train(no2~. , data, method = "rf", trControl = control, na.action = na.exclude)
fit.xg <- train(no2~. , data, method = "xgbLinear", trControl = control, na.action = na.exclude)
fit.gl <- train(no2~. , data, method = "glmnet", trControl = control, na.action = na.exclude)
results <- resamples(list(RF = fit.rf, XG = fit.xg, GL = fit.gl))
Y, por último, tendríamos que analizar qué modelo nos conviene tomar:
summary(results)
scales <- list(x = list(relation = "free"), y = list(relation = "free"))
bwplot(results, scales = scales)
Aquí hay más ejemplos del análisis para decidir entre los distintos algoritmos.